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ABSTRACT 


This study determines how to invest limited resources to increase the resilience of an 
infrastructure system against both non-deliberate and deliberate events. We propose an 
optimization model that seeks the best defensive investment for a weighted combination 
of deliberate and non-deliberate events. We formulate the general problem and conduct 
numerical analysis using a specific infrastructure system as a concrete example. We 
perform parametric analysis on the combined model in order to explore the way in which 
different solutions depend on the distributed weights and yield insight into the best 
investment decisions. 
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EXECUTIVE SUMMARY 


Owners and operators of critical infrastructure systems face a challenge when it comes to 
investing limited resources to maintain continuity of function in the presence of 
disruptive events. On the one hand, they need to protect against the most commonly 
occurring disruptions (e.g., failures or weather events). On the other hand, they need to 
protect against worst-case disruptions, as might be caused by an intelligent adversary 
who uses insider knowledge to deliberately interrupt system function. Often these two 
types of potential disruptions point to different parts of the system as “critical,” and in 
that case this tension can lead to a dilemma in terms of what to protect. 

In this study, we consider how to invest limited resources to increase the 
resilience of an infrastructure system against both non-deliberate and deliberate events. 
We represent the function of the infrastructure as a network flow problem and use an 
optimization model that seeks the best defensive investment for a weighted combination 
of deliberate and non-deliberate events. We conduct parametric analysis on the combined 
model in order to explore the way in which different solutions depend on the distributed 
weights and yield insight into the best investment decisions. From our analysis, we 
discover that it is possible to choose the best defense for that specific network by 
mapping the frontier between the worst-case and random disruptions. The results of that 
frontier gives us insight into (1) the minimum number of edges to protect, (2) the 
maximum defense for both worst-case and random disruptions, and (3) a moderate cost- 
effective defense for both worst-case and random disruptions 

Our investigation of a small network system provides proof-of-concept that the 
technique we propose can be effective, but there is more work to be done. In particular, 
our solutions are computationally intensive even for small systems, and additional 
research to reduce the required computational effort is needed before the technique can 
be used at large scale. 
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I. 


INTRODUCTION 


The United States (U.S.) has eritieal infrastrueture systems (e.g., eleetrie power, 
transportation, teleeommunioations, water) that need to be protected from disruptive 
events that are non-deliberate (e.g., natural disasters, accidents, failures) and deliberate 
(e.g., vandalism, sabotage, terrorism). Recently, these systems have been the targets 
terrorist attacks (September 11, 2001) and natural disasters (Hurricane Katrina in 2005 
and Hurricane Sandy in 2012) which caused billions of dollars in damage. In response to 
these and other disruptive events, the U.S. is spending billions of tax payer dollars 
annually (e.g., $50 billion in 2014) to strengthen infrastructure (The White House 2014). 
Of particular interest is how to allocate these funds in a manner that increases the 
resilience of our critical infrastructure systems. 

In the reliability engineering literature (e.g., Hwang et al. 1981 and Billinton and 
Allan 1992) and the risk literature (e.g., Bedford and Cook 2001), non-deliberate 
disruptions are commonly modeled as random events using probabilities. However, 
deliberate disruptions are more appropriately represented as worst-case events 
using game theory (e.g.. Brown et al. 2006). These two perspectives often point to 
different parts of an infrastructure system as “critical.” This leaves policy makers in a 
quandary about how to invest a limited defensive budget to make the infrastructure 
system more resilient. 

The goal of this study is to determine how to invest limited resources to increase 
the resilience of an infrastructure against both non-deliberate and deliberate events. 
This information will provide decision makers with insight into which assets to protect 
and how much emphasis to put toward preparing for natural disasters vice a 
terrorist attack. 

We propose an optimization model that seeks the best defensive investment 
for a weighted combination of deliberate and non-deliberate events. We formulate the 
general problem and conduct numerical analysis using a specific infrastructure system 
as a concrete example. We perform parametric analysis on the weights in the combined 
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model in order to explore the way in whieh different solutions depend on the distributed 
weights and yield insight into the best investment deeisions. 
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II. LITERATURE REVIEW 


Presidential Deeision Direetive 21 (PPD21), issued in February 2013, establishes 
national poliey on eritieal infrastrueture seeurity and resilienee (The White House 2013). 
Following the definition established previously in the USA PATRIOT Aet (Department 
of Justiee (DOJ) 2001), PPD21 defines eritieal infrastrueture to inelude “systems and 
assets, whether physieal or virtual, so vital to the United States that the ineapaeity or 
destruetion of sueh systems and assets would have a debilitating impaet on seeurity, 
national eeonomie seeurity, national publie health or safety, or any eombination of those 
matters.” 

PPD21 speeifieally states “It is the poliey of the United States to strengthen the 
seeurity and resilienee of its eritieal infrastrueture... to take proaetive steps to manage 
risk and strengthen the seeurity and resilienee of the Nation's eritieal infrastrueture, 
eonsidering all hazards that eould have a debilitating impaet on national seeurity, 
eeonomie stability, publie health and safety, or any eombination thereof These efforts 
shall seek to reduee vulnerabilities, minimize eonsequenees, identify and disrupt threats, 
and hasten response and reeovery efforts related to eritieal infrastrueture.” In its 
definition of “all hazards,” PPD21 refers to “a threat or an ineident, natural or manmade 
[that]... ineludes natural disasters, eyber ineidents, industrial aeeidents, pandemies, aets 
of terrorism, sabotage, and destruetive eriminal aetivity targeting eritieal infrastrueture.” 

In response to PPD2I, the Department of Homeland Seeurity (DHS) issued an 
updated National Infrastrueture Proteetion Plan (NIPP) in Deeember 2013. The 2013 
edition of the NIPP “eontinues to foeus on risk management as the foundation of eritieal 
infrastrueture seeurity and resilienee” (p. 4). 

The applieation of risk analysis to eritieal infrastrueture systems follows a long 
history in the applieation of probabilistie risk assessment (PRA) to eomplex engineering 
systems (see Bedford and Cook 2001 for an introduetion, or National Aeronauties and 
Spaee Administration (NASA) 2011 for a reeent applieation of PRA to spaee mission 
planning). Following the attaeks of 9/11, mueh of the work on risk to eritieal 
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infrastructure has focused on the application of PRA to terrorism risk (e.g., Pate-Comell 
and Guikema 2002, Garrick et al. 2004, Parnell, et al. 2005, Willis 2007). The 
Department of Homeland Security has promoted the use of probabilistic techniques when 
assessing the risk critical infrastructure systems (DHS 2013). 

A different line of research has focused on the ability of an intelligent adversary 
to deliberately disrupt system function. Specifically, there is now considerable work on 
the use of system interdiction models for understanding homeland security problems and 
the risk to critical infrastructure systems—these models are known as Attacker-Defender 
(AD) and Defender-Attacker-Defender (DAD) problems (Brown et ah, 2005, Brown et 
ah, 2006, Alderson et ah, 2014a, b). Applications of this model have been used in a 
variety of systems: fuel (Ileto 2011), cargo (De La Cruz 2011), the electric grid 
(Salmeron et al. 2004, 2009), and military networks (Burton 2013, Long 2013). 

In this thesis, we propose a model to defend our infrastructure against both non- 
deliberate hazards and deliberate threats. We use a linear combination of the two 
techniques previously described in order to map the efficient frontier of the tradeoffs 
between the consequences of a worst-case attack and expected consequence over a set of 
random events that will give decision makers the ability to choose the type of protection 
needed for his or her particular network. 

Previous work by Bier (2007), Zhuang and Bier (2007), and Hausken et al. (2009) 
has looked at the challenge of choosing what to protect, particularly when facing “all 
hazards” to include deliberate and non-deliberate events. However, these studies have 
abstracted the consequence of the event to a simple known function. As demonstrated by 
Alderson et al. (2013), in general it is difficult to approximate the consequences of losing 
one or more infrastructure components, or to pick the most vital components in a network 
system. An explicit objective of this research is to connect the operational models of 
infrastructure function and consequence assessment with the need to identify protective 
measures that work for both deliberate and non-deliberate disruptive events. 
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III. MODEL DEVELOPMENT 


In this chapter, we present a sequence of models that we use to perform our 
analysis. Throughout, we follow the concepts and notation in Alderson et al. (2014a). A 
speeific implementation of these models follows in Chapter IV. 

A. NORMAL OPERATION 

In this section, we provide an overview of the sequence of models to follow. Our 
starting point is an Operator Model of the form 

min/(w,i,y), (1) 

where w denotes the fixed design of the system, x denotes the operational setting (e.g., 
the status of each component of the system), and Y(w) denotes the set of feasible 
activities available to the system operator, given design w. 

The Operator Model represents the deeision making of an entity, called the 
operator, who ehooses the aetivities in the system. The Operator Model is prescriptive; its 
solution y* = argmin/(w,x,>’) represents the activities in the system that yield the 

_veT(w) 

minimum operating eost. 

The Operator Model (1) is a mathematical program that can be solved using 
various teehniques (see Alderson et al. 2014a for a detailed discussion). 

B. WORST-CASE DISRUPTIONS 

We are interested in understanding how well the system performs in the presenee 
of worst-case disruptions and/or most likely disruptions. For a fixed design w, we define 

P(w) = max min / {w,x,y ), (2) 

xeX yeY{w) 

where X is the set of feasible disruptive events, and x =argmax min f{w,x,y) is a 

worst-case disruptive event. Here, P(w) represents the performance of the system in the 

presence of the worst-case disruptive event. We refer to (2) as the Attacker Model; it 

builds on the previous formulation but has additional elements. 
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In general, P(w) can be evaluated three ways: 

1. by exhaustive enumeration of xgX , 

2. by taking the dual of the inner Operator Model to transform the max-min 
of the Attacker Model into a max-max problem (e.g., Brown et ah, 2006), 
and 

3. using Bender’s decomposition (e.g.. Benders, 1962). 

Given a set of possible designs W , our goal is to find a design wgW that yields 
“good” performance in the presence of disruptive events. 

If we are only concerned with defending against worst-case disruptions, then we 
solve a problem of the form: min P(w). This is simply a deterministic Defender- 

weW 

Attacker-Defender (DAD) problem. Let c^W denote the set of optimal designs that 
minimize P(w). That is, w* e w* = argminT’(w) . 

weW 

In principle, solving the Defender Model requires nothing more than enumerating 
every possible combination of defense and attack, then solving the corresponding 
Operator Model for each, and then finding the one that yields the lowest cost. This is the 
first half of our combined model which will determine the effect of the worst case. The 
second half of the model is the probabilistic model that optimizes defense against 
hazards. 

C. PROBABILISTIC DISRUPTION (MOST LIKELY) 

For a fixed design w, define 

Q(w) = E- min /(w, x, y ), (3) 

_V€7(w) 

where x is a random variable having a specified distribution, and E. denotes the 

expectation relative to x. Here, Q(w) represents the average performance of the system 
in the presence of a random disruptive event. 

In general, Q(w) can be evaluated two ways: 

1. via Monte Carlo sampling on x, and 
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2. using exhaustive enumeration if x is discrete to compute an exact 
expectation. 


If we are only concerned with defending against random disruptive events, then 
we solve a problem of the form: min Q{w). This is simply a stochastic optimization 

weW 

problem. Let W* eW denote the set of optimal designs that minimize Q(w). That is, 
w* G w* = argming(w) . 

weW 


D. DEFENDING AGAINST WORST-CASE AND MOST LIKELY 

In practice, we are concerned with both random disruptive events and worst-case 
disruptive events. To determine a good design, we must take into considerations the 
following two possibilities. 

1. Observe that if J¥^ ^ 0 , then we have an optimal solution to both 

problems and there is no issue. We simply choose w* e {W^ ■ 

2. However, in general we expect = 0 . More strongly, our 

experience is often that the worst-case disruptive event is very different 
from the average case disruptive event. This creates a dilemma in terms of 
where to invest defensive resources to make the system more resilient. 

We propose to study the combined problem as 

Z * (Z) = min[ZP(w) + (1 - Z)g(w)] where 0 < Z < 1. (4) 

wePF 

Note that when Z = 0, Eq.(4) corresponds to the stochastic average-case. In 
contrast, when Z = 1, Eq.(4) corresponds to the deterministic worst-case. 
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IV. CASE STUDY 


In this thesis, we use a small numerieal example to illustrate our teehnique. 
Following Alderson et al. (2014b), Figure 1 displays a notional (fuel distribution) 
infrastrueture system made of nodes and links. The system is a single eommodity 
transshipment model, designed to move fuel from supply loeations to demand loeations. 



Figure 1. Notional fuel distribution network (from Alderson et al. 2014b). 


In Figure I, the white cireles represent loeations with demand = 1 and the black 
circles represent locations with supply =10. All links are bidirectional and have per unit 
flow usage cost = 1. Each link has capacity =15 and the penalty for unsatisfied demand 
per node =10. Usage penalty for an interdicted arc = 20 per unit of flow (e.g., hiring a 
truck to drive the unit of flow to bypass the link). The objective is to minimize the sum of 
the flow cost for fuel deliveries and any penalties from unsatisfied demands. 

As shown in Figure 1, Nodes 3, 4, and 16 have two (parallel, redundant) 
connections to the rest of the network. This network has been built to be N-\ reliable, 
meaning that the loss of any single link does not disconnect any node. 


9 




A. OPERATOR MODEL 


We formulate the Operator Model for this system as follows. 

Indiees and Sets 

n EN nodes (alias i,j) 

[/j] EE undireeted edge between nodes i and j, i<j 

{i,j) EA direeted are from node i to node j 

[i,j] EE (/</•) A {(ij) EA A ij,i) EA)) 
dED defenses, Z)={0,1}; 0 represents undefended, 1 represents 

defended (and invulnerable) 

Data [units] 

Cij per unit eost of traversing are {i,j) EA [dollars/barrel] 

Uij upper bound on total (undireeted) flow on edge [iJ] EE 

[barrels] 


Xij 1 if edge [i,j] EE damaged, 0 otherwise [binary] 

q‘^. per unit penalty eost of traversing are (i,j) EA if damaged 
[dollars/barrel] 

bn fuel supply at node n EN [barrels] 

(-demand for bn <0 ) 

Vn per unit penalty eost for demand shortfall n EN 

[dollars/barrel] 

w‘‘j 1 if defense option d is implemented on edge [i,j] EE, 

0 otherwise [binary] 
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Decision Variables [units] 


flow on arc {i,j) EA for defense d [barrels] 
Sn fuel shortfall at node n EN [barrels] 


Formulation; 


min 

Y,S 

Z Z + (s- + ) Y' 

cIeD [iJ]eE 

] + Z^A 

neN 

(DO) 

s.t. 

z z k-z z Y,:-s„<b, 

deD j:(nJ)EA deD i:(i,n)eA 

yne N 

[^„] (Dl) 


U Ji u U 

V[z, A,d ^ D 

[^4] (D2) 


o 

Al 

V[z, j]^ A,d G D 

(D3) 



yne N 

(D4) 


Discussion 

The objective function (DO) combines the total flow cost and the total penalty 
cost. Constraints (Dl) enforce balance of flow at each node. Stipulations (D2) and (D3) 
ensure bounds on decision variables. This formulation implements cost-based 
interdiction—that is, damage to an arc makes it extremely expensive but not infeasible— 
which makes the problem easier to solve computationally. In the above example, we have 
b8=bio=\Q and h„=-l otherwise. In addition, we set Cij=\, q\ =Q,y{i,j) ^ A 

u.j =100,V(/,7)e = 10, V(z, 7 ) e , Uij=\5 for all [zj] EE, and v„ =10 for all nEN. 

There are no arcs defended in the base case; therefore, =0 andw° =1 for all [zj] EY. 

Solving this problem, we obtain a minimum-cost solution consisting of flows that 
meet all the demands (see Figure 2), with the flows labeled in blue. Our findings are 
consistent to the base solution findings in Alderson et al. (2014b). We are concerned with 
the ability of the system to perform in the presence of either random (probabilistic) or 
deliberate (worst-case) disruptive events. 


11 




Figure 2. Operator Model base solution 


B. ATTACKER MODEL 

For this notional fuel system, we model the attacker’s problem as follows. 


Additional Data[units] 

Tij “cost” to break edge [zj] E E [cardinality] 

attack_budget budget constraint on the number of simultaneous attacks 

[cardinality] 


Additional Decision Variables [units] 

Xij 1 if attacker breaks edge [z j] EE, 0 otherwise [binary] 


Formulation 
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P{w) = max min [(c.. + + (c 

rfeZ) [1,./]££■ 

neN 

(ADO) 

s.t. (D1),(D2),(D3)(D4) 



^ rijXij < attack _ budget 


(ADI) 


y[i,j] eE 

(AD2) 


The objective function (ADO) is the same as that for the Operator Model (DO), 
except that parameters Xij have been replaced by decision variables Xy. Constraint (ADI) 
limits the number of simultaneous attacks, and the cost to attack each arc can be different. 
Stimpulations (AD2) require that attacks are binary. We note that qtj = 0 implies that arc 
{i,j) is effectively invulnerable, because attacking it does not increase the flow cost for the 
operator. 

In the above example, we model parallel arcs as costing twice as much to attack. 
That is, we have r 2,3 = r 4_8 = r 12,16 = 2 and all other = 1. 

The solution for the Attacker Model depends on the attacker budget (the number 
of attacks that the attacker can afford). The attacker seeks to maximize the cost to the 
operator. In Figure 3, an attack on edge [10,13] represents the worst case attack to the 
operator when the attacker can afford to target only a single edge. The operator reroutes 
his flow as shown in blue and the overall operator cost increases from 25 to 33. An 
attacker having the ability to target a large number of edges could render our network 
practically useless as shown in Figure 4. Here, the attacker has the budget to target five 
edges and focuses on the arcs around the supply nodes to cripple flow to the entire 
network with exception to node 4. The attacker picks the best place to attack depending 
on how many attacks allotted by the attacking budget. These results are also consistent 
with Alderson et al. (2014b). 
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Figure 3. Operator Model solution with 1 attack 



Figure 4. Operator Model solution with 5 attacks 


C. DEFENDING AGAINST A WORST-CASE ATTACK 

In this thesis, we restrict defensive actions to the “hardening” of edges, such that 
targeting an edge does not increase its usage cost. We formulate the following model. 

Additional Data[units] 

defense_budget budget constraint on the number of defenses [cardinality] 
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Additional Decision Variables [units] 



1 if defense option d is implemented on edge [i,j] EE, 
0 otherwise [binary] 


Formulation 


mmmaxminX S +iCj, 

’ rfe£)[i,7]e£ 


(DADO) 

s.t. (D1)(AD1)(AD2)(D3)(D4) 




V[z,7] G E,d e D 

{DADl) 

^ Wl < defense_budget 


(DADl) 

II 

V[i,j]GE 

{DAD3) 

W‘ E {0,1} 

V[z,7] G E,d gD 

{DAD4) 


The objective (DADO) includes the cost of flow over existing (and possibly 
damaged or protected) edges and penalties for unmet demand. Contraints (DADl) 
ensures that flow on edge [i,j] does not exceed flow capacity Wy. The constraint (DAD2) 
requires that the cost of all defenses fall within the existing defense budget using the 
binary variable . Constraints (DADS) ensure that flow on each edge \i,j~\ is on either 

the protected or unprotected edge and not both. Here, the pair Wlj=\ and W°=Q denotes a 

defended edge and vice versa for undefended. Stipulations (DAD4) require binary 
defense decisions. 

In order to be consistent with Alderson et al. (2014b) we will give examples of 
defense with an attacker budget of three attacks as shown in Figure 5. For the purposes of 
this study we only focus on hardening edges. 
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Figure 5. A. AD Model: worst-case attack involving three edges. In the absence 
of any defenses, edges [2,7], [10,13], and [11,15] are targeted. B. DAD 
Model: the corresponding worst-case attack with one edge protected. 

Edge [10,13] is defended, and edges [6,10], [7,8], and [9,13] are 
targeted. 

The defender model protects edge [10, 13], note that it is the same arc the attacker 
chose as the one arc attack in Figure 3. The attacker picks another set of arcs to attack in 
order to maximize cost for the operator that does not include the defended arc. With edge 
[10,13] defended and in the presence of three targeted edges, the operator is able to 
satisfy demand at seven nodes, instead of only six with no arcs defended. 

The defender can limit the attacker by strengthening more arcs at shown in Figure 
6. However, we cannot defend all arcs because of budget constraints. Therefore, 
decisions about what arcs to protect must be made within those constraints to minimize 
cost to the operator. 
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Figure 6. Worst-case disruption with 3 edges targeted and 5 edges protected. 

D. EXPECTATION MODEL 

We use the symbol x to represent a random event in place of x, which denotes a 
set of edges that have been damaged. We calculate the expected cost due to random 
disruptions, and will not focus on a single, worst-case event. For example, 
S = |x| ^x, < 2| represents the set of all events with two or fewer li nk s broken. For any 

element xeS, let be the probability that x occurs, i.e., {0</f<l VxeS , 
^ if = 1}. For simplicity, we assume the probability that an edge will go down is 

xeE 

independent between edges, giving it a binomial distribution with x becoming a random 
variable. We can think of S as a list of the enumerated (all) possible attacks with 
associated probabilities. 

Additional Sets 

e S set of (enumerated) outcomes 

Additional data [units] 
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c {0,1}'^' set of binary vectors, one element for each edge, where x| =1 
if edge [z, 7 ]e£'is damaged, and =0 if the edge remains intact, in 
outcome e S . 


Formulation 

Q{w) = E^[f[w,^)) , where 

F (w,= mm X Z [(% + + (S.- + ] + X ^nS„ 

d&D[i,j'\&E n^N 

s.t. (D3)(D4) 

Z Z K-H Z Yf„-Sn<bn \/n^N 

d^D j\{n,j)^A d&D 

^ij + V[z,y] G J G Z) 


(EO) 


[^f] (El) 
[F,f ] (E2) 


Whereas the AD model yields both the worst-case attack and associated best- 
response in the objective value, the result of the ED model is solely the objective value. 

Defending against the average case disruption means we must take into account 
the arcs most likely to fail and either strengthen or build around those arcs in each x . 

E. DEFENDING AGAINST THE EXPECTED CASE 

Just like in the DAD model, we are defending edges in order to minimize the 
system cost in the presence of disruption by carefully choosing which edges to protect. 

The possible defenses will be the same as in the DAD model, just like in Figure 7. 
However, the result will be the objective value comprised of the products from this 
design and all possible failures in E creating the defender-expectation-defender (DED) 
case. 
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Figure 7. Example Defense in which edges [2,7], [7,8], [9,13], [10,11], and 
[10,13] are defended. 


Formulation 

In order to determine the best defense for the expected value problem, we use 
Benders decomposition (Benders, 1962), an iterative algorithm in which we use the index 
k to denote an iteration, and K to denote the current iteration of the algorithm. The 
subproblem for any fixed defense w and any specific outcome e S is just 

and can be solved as a linear program with associated duals as given in (El) 

and (E2). Each of these dual solutions provides a linear lower bound on the objective 
function value of as a function of the defense: 

n&N [iJ]&Ed^D 

The weighted sum of these linear functions is thus a lower bound on the expected 
cost as a function of defense, so we then define: 

and 

^gH 

to be the weighted combinations of the optimal dual solutions to each of the subproblems 


19 

























solved in iteration k, yielding a single cut for the set of constraints (SMI) in 
the master problem. We therefore have the following stochastic master problem: 


min 


(SMO) 

S-1- + Z 

n&N [iJ^^Ed^D 

k = l,...,K 

(SMI) 

d^D 

V[/,7]eE 

(SM2) 

^ W^j < defense_budget 


(SMS) 


The objective (SMO) represents the expected cost of operating in the face of 
random events. We minimize this objective through each iteration k as shown in (SMI). 
Each edge [fy] has two possibilities, denoted with d, and only one can be used at a time 
shown in constraint (SM2). The constraint (SMS) ensures that the defenses used are 
limited to defense_budget. We describe the procedure for solving this problem in the 
context of the combined model, which now follows. 

F. COMBINED MODEL 

If we define to be the weight given to the stochastic (i.e., expected-value) 

objective, where 0 < < 1, and wt^^ = \- to be the remaining weight placed on 

the deterministic (i.e., worst-case) objective, then we can formulate the following 
combined master problem: 
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min 

z,w 


(CMO) 


s.t. 


^^stoch'^stoch ^^def^det 


(SM1),(SM2),(SM3) 


2de, Z [(Cy 


d^D [i,y]e£ 


nE.N 

k = \,...,K 

(CM2) 

Z Z k*-z z c 

deD j^nJ^eA JgZ) 

-Sl<b^ yneN,k = l,...,K 

(CMS) 

ydk Y‘i'^ < u W!^ 

u Jt ij 

k/[i,j] e E,d e D,k = 1,.. 

.,K (CM4) 

o 

Al 

V(/,7) e A,d e D,k = \,. 

..,K (CM5) 

Co 

IV 

O 

\/n^N,k = \,...,K 

(CM6) 


Here the defense decisions, W, are made before any realization of an attack or a random 
failure, and the evaluation of the defense chosen at iteration A: is a convex combination of 
the evaluation of each of the two separate subproblems. 

At each iteration K of Benders decomposition, our algorithm solves each of the 
two subproblems (stochastic and deterministic) to optimality for a fixed defense, . 
The constraints from the stochastic master problem transfer directly to this combined 
master problem. The deterministic subproblem yields a single worst-case attack, , and 
this provides a lower bound on the cost of the optimal response for that particular attack, 
although only in terms of the flow variables in response to that attack, and this 
contributes a new cut to the set of constraints (CM2) in the master problem. We therefore 
must also create a block of nonnegative variables, , to model that response, and two 
additional blocks of constraints to ensure they (a) conform to balance of flow (CMS), and 
thus act like admissible flows in our network, and (b) correspond appropriately to any 
defense, W (CM4). The algorithm then solves the master problem for iteration K, and 
generates a new defense, . After the subproblems have been solved, we have an 
accurate evaluation of the combined cost of the current defense plan, . If this cost is 
lower than any seen so far, we record as a new incumbent solution and take its cost 
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as the new upper bound. Eaeh solve of the master problem provides a lower bound on the 
optimal solution, and these values are monotonieally nondeereasing. If, at any point, the 
relative differenee between the best upper bound derived from the subproblems and the 
best lower bound from the master problem is less than a predetermined optimality gap, a, 
the algorithm terminates and reports the eurrent ineumbent as an ^■-optimal solution. This 
eombined model uses the duals from the stoehastie model and eombines them with the 
worst ease model to produee outs. 
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V. RESULTS AND ANALYSIS 


In this chapter, we explore the differenees between the “worst-ease” and “most 
likely” disruptions for the notional fuel distribution network in Figure 1. We use the 
eombined model to perform our analysis, where in equation (CMO) we have eonvex 
weights wtdet = ot and wtstoch = 1 - a. 


A. ANALYSIS OF FIVE SCENARIOS 

We ran the eombined model for five seenarios, eaeh deseribed in the subseetions 

below. 


I. Two Attacks and One Defense 

Figure 8 displays the results of the eombined model for our notional fuel system 
in the ease of two attaeks and one defense. For different values of a, we solve the 
eombined model to within an optimality toleranee of 0.001. On a eomputer with a 
dual-eore 1.2 GHz proeessor, the solution time for all seenarios involving two 
attaeks was approximately 2 hours for all eleven eases (a=0.0 through a=1.0). For eaeh 
value of a. Figure 8 shows the value for and along with their eombination 

weighted by a . 

When defending solely against the worst-ease attaek (a=l), the best defense is 
[2,7] and the worst two-link attaek is [7,8] and [8,12]. When defending solely against 
random events (a=0), the best defense is also [2,7]. Note that the value Z^^^ is 
signifieantly higher than meaning that the worst-ease disruption is signifieantly 

worse than the average ease. Beeause the worst-ease solution does not vary with a, the 
values for Z^^^^^and are insensitive to ehanges in a. As a result, the eombined 

objeetive value is not meaningful on its own, rather it represents simply the weighted 
eombination of these two eonstant values. 
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Weight (a) 


•Zdet 

•Zstoch 

■Objective_(CMO) 


Figure 8. The objective values for two attacks and one defense versus the weight 
for given weight a. 


2. Two Attacks and Two Defenses 

Just as in the previous case, the worst-case solution does not vary with a, the 
values for and are insensitive to changes in a. When defending solely against 

the worst-case attack (a=l), the best defense is [2,7] and [7,8]; while, the worst two-link 
attack is [1,2] and [9,13]. When defending solely against random events (a=0), the best 
defense is also [2,7] and [7,8]. 

Table 1 summarizes the results from the scenarios involving no more than two 
broken links. 


Table 1. Optimal links to defend for different defensive budgets for cases 
when there are no more than two broken links. For this small 
network, we observe that the single best link to protect is [2,7], and 
when defending two links the best two are [2,7] and [7,8]. 



Stoch 

Combined 

Det 

Def 

(a=0) 

(a = 5) 

(a=l) 

1 

[2,7] 

[2,7] 

[2,7] 

2 

[2,7][7,8] 

[2,7][7,8] 

[2,7][7,8] 
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3. 


Three Attacks and One Defense 


The results of the combined model for our notional fuel system in the case of 
three attacks and one defense has a shape similar to the two-attack case. For this scenario, 
we solve the combined model for different values of a to within an optimality tolerance 
of 0.1. On a computer with a dual-core 1.2 GHz processor, the solution time for was 
approximately 6 hours for all eleven cases (a=0.0 through a=1.0). Table 2 displays the 
results of this computation. 


Table 2. Combined model results for three attacks and one defense 


a 

edge to defend 

worst attaek 

objective (CMO) 

obj Zdet 

obj Zstoch 

0.0 (Stoch) 

12,71 

[9,131[6,101[7,81 

29.47955282 

80 

29.47955282 

0.1 

12,71 

[9,13116,10117,81 

34.53159754 

80 

29.47955282 

0.2 

12,71 

19,13116,10117,81 

39.58364226 

80 

29.47955282 

0.3 

12,71 

19,13116,10117,81 

44.63568697 

80 

29.47955282 

0.4 

12,71 

19,13116,10117,81 

49.68773169 

80 

29.47955282 

0.5 

12,71 

19,13116,10117,81 

54.73977641 

80 

29.47955282 

0.6 

12,71 

19,13116,10117,81 

59.79182113 

80 

29.47955282 

0.7 

12,71 

19,13116,10117,81 

64.84386585 

80 

29.47955282 

0.8 

12,71 

19,13116,10117,81 

69.89591056 

80 

29.47955282 

0.9 

113,101 

19,13116,10117,81 

75.02040504 

80 

30.2040504 

1.0 (Del) 

113,101 

[9,13][6,101[7,81 

80 

80 

30.2040504 


When defending solely against the worst-case attack (a=l), the best defense is 
[13,10] and the worst three-link attack is [7,8], [9,13] and [6,10] as shown in Table 2. 
While defending solely against random events (a=0), the best defense is [2,7]. There is a 
difference in the solution for the worst-case and the solution for random events; however, 
the solution doesn’t change until a > 0.9 toward the worst-case. The difference in cost to 
the operator between protecting edge [2,7] and [13,10] is one. Similar to the two attack 
scenarios, the value is significantly higher than There is more value in 

protecting [2,7] for the stochastic case than [13,10], but the values for Z^^^_,^are similar for 
the two cases. Overall, the values for and are insensitive to changes in a, but 
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are far from each other. Again, the combined objective value represents simply the 
weighted combination of these two constant values. 


4. Three Attacks and Two Defenses 

We similarly solve for the results of the combined model for our notional fuel 
system in the case of three attacks and two defenses. On a computer with a dual-core 1.2 
GHz processor, the solution time for was approximately 10 hours for all eleven cases 
(a=0.0 through a=1.0). Table 3 displays the results of this computation. 


Table 3. Combined model results for three attacks and two defenses 


a 

edge to defend 

worst attaek 

objeetive (CMO) 

obj Zdet 

obj Zstoeh 

0.0 (Stoch) 

12,7117,81 

11,21113,101111,151 

28.03713818 

72 

28.03713818 

0.1 

12,7117,81 

11,21113,101111,151 

32.43342436 

72 

28.03713818 

0.2 

113,10117,81 

113,14118,121110,111 

36.47357923 

67 

28.84197404 

0.3 

113,10117,81 

113,14118,121110,111 

40.28938182 

67 

28.84197404 

0.4 

113,10117,81 

113,14118,121110,111 

44.10518442 

67 

28.84197404 

0.5 

113,10117,81 

113,14118,121110,111 

47.92098702 

67 

28.84197404 

0.6 

113,10117,81 

113,14118,121110,111 

51.73678961 

67 

28.84197404 

0.7 

113,10117,81 

113,14118,121110,111 

55.55259221 

67 

28.84197404 

0.8 

113,10117,81 

113,14118,121110,111 

59.36839481 

67 

28.84197404 

0.9 

113,10117,81 

113,14118,121110,111 

63.1841974 

67 

28.84197404 

1.0 (Del) 

113,10117,81 

113,14118,121110,111 

67 

67 

28.84197404 


The difference between defending solely against random events (a=0) and the 
worst-case attack (a=l) is the defending of edge [2,7] vice [13,10]. Although [13,10] and 
[2,7] yield the two highest operating costs for attacking a single edge, it is more 
important to protect only one of them along with [7,8] due the construction of the 
network as shown in Figure 1. Unlike the previous defense for three attacks, the solution 
turns to the worst-case when a > 0.2, meaning, when coupled with the defense of [7,8] it 
is more beneficial for both models to protect [13,10] vice [2,7]. 

5. Three Attacks and Three Defenses 

We similarly solve for the results of the combined model for our notional fuel 
system in the case of three attacks and three defenses. On a computer with a dual-core 1.2 
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GHz processor, the solution time for was approximately 13 hours for all eleven eases 
(a=0.0 through a=1.0). Table 4 displays the results of this eomputation. 


Table 4. Combined model results for three attaeks and three defenses 


a 

edge to defend 

worst attack 

objective (CMO) 

obj Zdet 

obj Zstoch 

0.0 (Stoch) 

[2,7117,8119,131 

11,21113,101111,151 

27.45686501 

72 

27.45686501 

0.1 

12,7117,81113,101 

18,121113,141110,111 

31.26976353 

67 

27.29973726 

0.2 

12,7117,81113,101 

18,121113,141110,111 

35.23978981 

67 

27.29973726 

0.3 

113,10117,8118,121 

12,7119,131110,111 

39.12518556 

64 

28.4645508 

0.4 

113,10117,8118,121 

12,7119,131110,111 

42.67873048 

64 

28.4645508 

0.5 

113,10117,8118,121 

12,7119,131110,111 

46.2322754 

64 

28.4645508 

0.6 

113,10117,8118,121 

12,7119,131110,111 

49.78582032 

64 

28.4645508 

0.7 

113,10117,8118,121 

12,7119,131110,111 

53.33936524 

64 

28.4645508 

0.8 

113,10117,8118,121 

12,7119,131110,111 

56.89291016 

64 

28.4645508 

0.9 

113,10117,8118,121 

12,7119,131110,111 

60.44645508 

64 

28.4645508 

1.0 (Det) 

113,10117,81110,111 

12,7119,13118,121 

64 

64 

28.68990798 


The three attaek and three defense solution for the eombined model isn’t as 
intuitive as the previous two seenarios. When0.3 <a< 0.9 , the optimal solution defends 
edges [8,12], [13,10], and [7,8]. Showing that moderate proteetion against both random 
and worst-ease disruptions is separate from the solutions for either event as shown in 
Table 2. 

Figure 9 shows that as the values for the separate models ehange, the objeetive 
funetion still maintains the linear eharaeteristie as expeeted; however the operator eosts 
for Zjgjgoes down for eaeh ehange in defense while goes up slightly for the 

inerease in a. The worst ease and expeetation models both ehose the same edges to 
defend for the two attaek seenario as shown in the Table 1. This result makes it easy for a 
deeision maker to ehoose whieh items to defend for the best possible proteetion. The 
problem for deeision makers is shown in the ease of three attaeks. 
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Weight (a) 


•Z_det 

•Z_stoch 

•Objective_(CMO) 


Figure 9. The objective values for three attacks and three defenses versus the 
weight for given weight a. 


B. DISCUSSION 

In the case of random disruptions, with a limited number of defenses, we might 
not be able to protect our system against the (potentially enormous) set of disruptions 
near the “center” of the distribution. However, if there are only a few disruptions that 
cause significant damage (i.e., well above the mean) then we can apply our limited 
resources to reduce their impact and therefore have a noticeable effect on the expected 
cost over the entire distribution of random disruptions. 

Figure 10 displays the overall value of having one, two, or three defenses, when 
protecting against three attacks. There are diminishing returns for the third defense used. 
To save money the decision maker should use a two defense approach vice three in this 
case. 


28 












Weight (a) 


)( '0 defense' 



defense' 

defense' 


•'3 


defense' 


Figure 10. Multiple defense options against the three attacks scenario; Objective 
value for given defense with weight a in equation (CMO), where wtdet 
= a and wtstoch = 1 - a. 
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VI. CONCLUSIONS 


In this thesis, we have presented a model for defensive investment that eombines 
the threat of a worst-ease disruption with the risk of a random disruption. Our 
eombination model has the potential to help deeision makers come up with a good 
defense that fits their particular network. From our analysis we discover that it is possible 
to choose the best defense for that given network by mapping the frontier between the 
worst-case and random disruptions. 

Our investigation of a small network system provides proof-of-concept that the 
technique we propose can be effective, but there is more work to be done. We discovered 
that the stochastic model runs longer than the combined model. Interestingly, this 
suggests that knowing how to defend against the worst-case potentially provides insight 
into how to defend against random events. However, even for the small example here, the 
program takes a long period of time to solve the combined model for three attacks and 
two or more defenses. Therefore, additional research to reduce the required 
computational effort is needed before the technique can be used at large scale. 
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APPENDIX 


Our computational models and solution were obtained with the Python 
programming language and the Cooper-Pyomo optimization module using CPLEX as the 
solver. The program is solved using five modules. The combo_main. py module 
contains the outer loop for multiple runs of the model and is the driver that branehes into 
the other modules. We have a separate module, phat_code .py, to create the 
distribution for the probability of damaging one or more arcs. The network is built using 
the networkx module, in netx_graph . py, by building the list of edges and nodes. The 
file combo_model. py creates the pyomo models for the master and sub problems used 
to create cuts and implement benders. The file pyomo_model .py creates and solves 
the minimum cost flow model. The code for the modules are listed below: 


#Begin combo_main.py 

import netx_graph 

import pyo_model 

import combo_model 

import phat_code 

import matplotlib.pyplot as pit 

#open file 

file= open ( 'ComboResults.txt' , 'w') 

file.write('Weight (Worst), Arcs Defended, Objective value. Worst attaclc,Worst 
value. Random valueXn') 

test = [0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,11 
for weight in test: 

if _name_ == '_main_' : 


#establish nodes, edges, start_nodes, and double_edges 

#This will eventually occur in a subroutine that parses input files 

nodes = ['1', '2', '3', '4', '5', '6', '7', '9', 'll', '12', '13', '14', '15', '16'] 

edges= 

[('!', '2'), ('2', '7'),('!', '5'), ('5', '9'), ('9', '13'), ('13', '10'), ('13', '14'), ('1 
4', '15'), ('6', '10'), ('6', '7'), ('7', '8'), ('8', '12'), ('10', 'll'), ('11', '12'), ('ll 
', '15') ] 

start_nodes = ( '10 ' , ' 8 ' ) 

double_edges =[('2', '3'), ( ' 8 ' , ' 4 ' ) , ( ' 12 ' , ' 16 ' )] 

#build graph 

Grph, vuln, double_nodes = netx_graph.newmodel(nodes, start_nodes, 
edges, double_edges) 
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#establish single-edge probabilities and initialize xhat 

phat, xhat = phat_code.phat_edges(Grph, vuln) 

defenses = ['O', '1'] 

w_hat = {} 

for (i,j) in vuln: 

w_hat[i,j,'0'] = 1 
w_hat[i,j,'1'] = 0 


#Create dictionaries for node supplies and vulnerable edge capacities 
supply = {} 

for item in Grph.nodes () : 
if item in start_nodes: 

supply[item] = 10 
elif item in double_nodes: 

supply[item] = 0 
else: 

supply[item] = -1 
capacity = {} 
for (i,j) in vuln: 

capacity[i,j] = 10 

#generate probability map for all events with two or fewer failures 
P = phat_code.create_Prob(Grph,vuln,phat,3) 
print 'Number of events: ' + str(len(P)) 
max_iterations = 50 

nodepi = {} #dictionary of node duals, by iteration 
edgemu = {} #dictionary of arc capacity duals, by iteration 
xhats = {} #dictionary of worst-case attacks, by iteration 


LB=0 # can set this based on first MP solve... 

UB=10000 # should come up with a better way of handling this 

master_tol =0.10 # relative gap tolerance for terminating benders 

det_wt = weight 
stoch_wt = 1-weight 

for iteration in range(max_iterations): 
expected_cost = 0 
cumul_p = 0 

for (scen_num,x_tilde) in enumerate(P.keys ()) : 

xhat = phat_code.generate_xhat_from_x_tilde(vuln,x_tilde,xhat) 
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model, instance 

PYO_model.build_model (Grph, xhat,w_hat,vuln,defenses,supply,capacity) 

results = pyo_model.solve_model(instance) 

instance.load(resuits) 

objective_value = instance.FlowCost () 

if scen_num == 0: 

worst_scen = x_tilde 

worst_obj = objective_value 

for node in Grph.nodes () : 

dualval 

instance.dual.getValue(instance.FlowBalance[node]) 

if abs(dualval) > le-08: 

nodepi[node,iteration] = P[x_tilde]*dualval 
else: 


nodepi[node,iteration] = 0.0 
for (i,j) in vuln: 

for defense in defenses: 
dualval= 

instance.dual.getValue(instance.EdgeCapacity[i,j,defense]) 

if abs(dualval) > le-08: 

edgemu[i,j,defense,iteration]= 


P[x_tilde]*dualval 


else: 


else: 


edgemu[i, j,defense,iteration] = 0.0 


if worst_obj < objective_value: 
worst_obj = objective_value 
worst_scen = x_tilde 
for node in Grph.nodes () : 
dualval= 

instance.dual.getValue(instance.FlowBalance[node]) 

if abs(dualval) > le-08: 

nodepi[node,iteration] += P[x_tilde]*dualval 
for (i,j) in vuln: 

for defense in defenses: 
dualval= 

instance.dual.getValue(instance.EdgeCapacity[i, j, defense]) 

if abs(dualval) > le-08: 

edgemu[i,j,defense, iteration] += P[x_tilde]*dualval 


expected_cost += P[x_tilde] * objective_value 
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cumul_p += P[x_tilde] 

#provide intermittent updates 
#if scen_num % 25 == 0: 

# print str(scen_num) + ': ' + str(objective_value) + ' 

approx, expected cost: ' + str(expected_cost/cumul_p) 

#if len(x_tilde)==0: 

# print '***' + str(scen_num) + ' : ' + str(objective_value) + 

' approx, expected cost: ' + str(expected_cost/cumul_p) 

print 'Subproblem Iteration:' + str (iteration) 

print ' Expected Cost = ' + str(expected_cost) 

print ' Worst-Case Cost = ' + str(worst_obj) 

for (i,j) in worst_scen: 

print ' Attack: (' + str(i) + ', ' + str(j) + ') ' 

print ' Weighted Cost = ' + str(stoch_wt*expected_cost + 

det_wt *worst_ob j) 

#print edgemu 

#print nodepi 

if iteration==0: 

incumbent = w_hat 

UB = stoch_wt*expected_cost + det_wt*worst_obj 
inc_iteration = 0 

print '***new incumbent: (UB) ' + str(UB) 
x_star=worst_scen 
x_star_obj= worst_obj 
x_tilde_obj=expected_cost 
else: 

if UB > stoch_wt*expected_cost + det_wt*worst_obj: 

UB = stoch_wt*expected_cost + det_wt*worst_obj 

incumbent=w_hat 

inc_iteration = iteration 

print '***new incumbent: (UB) ' + str(UB) 
x_star=worst_scen 
x_star_obj= worst_obj 
x_tilde_obj=expected_cost 
if UB-LB<master_tol*LB: 
break 

for (i,j) in vuln: 

if (i,j) in worst_scen: 

xhats[i,j,iteration] = 1 
else: 
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xhats[i,j,iteration] 


0 


model, instance = combo_model.BuildComboMaster(Grph, vuln 
defenses, 3, iteration+1, det_wt, stoch_wt, xhats, nodepi, edgemu, supply 
capacity) 

results = combo_model.SolveModel(instance) 
instance.load(resuits) 

print 'Master Iteration: ' + str(iteration) 

w_hat = {} 

for (i,j) in vuln: 

for d in defenses: 

if instance.W[i,j,d].value > le-08: 
w_hat[i,j, d] = 1 
if d=='1': 


+ str(w_hat[i,j,d]) 


print 'Defend (' + str(i) + ',' + str(j) + ') W= 


else: 


w_hat[i,j, d] = 0 

if LB < instance.MasterObjective () : 

LB = instance.MasterObjective() 

print 'Master objective: (LB) ' + str(instance.MasterObjective()) 
if UB-LB<master_tol*LB: 
break 

print 

if UB-LB<master_tol*LB: 

print 'Optimality Gap Reached at iteration=' +str (iteration) 
else: 

print 'Iteration limit reached at iteration=' +str (iteration) 
print 'Incumbent Defense found at iteration=' +str(inc_iteration) 
for (i,j) in vuln: 

if incumbent[i,j,'l'] > le-08: 

print ' Defend (' + str(i) + ',' + str(j) + ')' 
file.write]'('+str(i) + ',' + str(j) + ')' ) 

print 'Expected Cost: ' + str(UB) 
print 'Lower Bound: ' + str (LB) 

if LB>0: 

print 'Optimality gap: ' + str((UB-LB)/LB) 

file.write(str(weight)+','+str(UB)+','+str(x_star)+','+str(x_star_obj)+','+str( 
x_tilde_obj)+'\n') 

file.close () 

#End combo_main.pY 
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#Begin netx_graph.py 

#Build the graph using networkx 

import networkx as nx 

def twolines(a,b,grph, vulnerable_edges, double_nodes): 

#create the equivalent of two parallel edges between nodes a and b by 
adding 

#four artificial nodes and twelve directed arcs 
suffix_l = + a + + b 

suffix_2 = + b + + a 

grph.add_node(a + suffix_l, node_size = 99,node_color='y') 

grph.add_node(b + suffix_l, color='y') 

grph.add_node(a + suffix_2, color='y') 

grph.add_node(b + suffix_2, color='y') 

grph.add_edge(a, a + suffix_l, cost =0) 

grph.add_edge(a, a + suffix_2, cost =0) 

grph.add_edge(a + suffix_l, a, cost =0) 

grph.add_edge(a + suffix_2, a, cost =0) 

grph.add_edge(b, b + suffix_l, cost=0) 

grph.add_edge(b, b + suffix_2, cost=0) 

grph.add_edge(b + suffix_l, b, cost=0) 

grph.add_edge(b + suffix_2, b, cost=0) 

grph.add_edge(a + suffix_l, b + suffix_l, cost = 1) 

grph.add_edge(a + suffix_2, b + suffix_2, cost = 1) 

grph.add_edge(b + suffix_l, a+suffix_l, cost = 1) 

grph.add_edge(b + suffix_2, a+suffix_2, cost = 1) 

vulnerable_edges.append( (a+suffix_l,b+suffix_l) ) 

vulnerable_edges.append( (a+suffix_2,b+suffix_2) ) 

double_nodes.append(a+suffix_l) 

double_nodes.append(a+suffix_2) 

double_nodes.append(b+suffix_l) 

double_nodes.append(b+suffix_2) 

return 

def doubs(a, b, grph, vulnerable_edges): 

#create two directed arcs to represent an undirected edge between nodes a 
and b 

grph.add_edge(a, b, cost=l) 
grph.add_edge(b, a, cost=l) 
vulnerable_edges.append( (a,b) ) 

return 

tCreate the model 

def newmodel(nodes, start_nodes, edges, double_edges): 
vulnerable_edges = [] 
double_nodes = [] 
g = nx.DiGraph() 

#node list 
for i in nodes: 

if i in start_nodes: 

g.add_node(i, node_color ='b') 
else: 

g.add_node(i, node_color='r') 

#edge list 

for i,j in edges: 

doubs(i,j,g,vulnerable_edges) 
for i,j in double_edges: 

twolines(i,j,g,vulnerable_edges, double_nodes) 
return g, vulnerable_edges, double_nodes 
#End netx_graph.py 
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#Begin phat_code.py 
timport random functions 
import random 

from itertools import combinations 
import math 

#Take a list of vulnerable edges and assign probabilities to them 
def phat_edges(graph,vulnerable_edges): 
random.seed(0) 
phat = {} 
xhat = { } 

#Initialize probablity of arcs with random and initialize xhat =0 
#for tups in graph.edges(): 
for tups in vulnerable_edges: 
temp = random.random() 
phat[tups]= temp*.3 + .01 
xhat [tups]= 0 
return phat, xhat 

def create_Prob(graph,vulnerable_edges,phat,max_losses=None): 

"""Build dictionary of events and their associated probabilities. 
Assumes elements (vulnerable_edges) fail independently with 
probabilities specified by the dictionary phat. If 
max_losses is specified, computes conditional probabilities 
for events whose cardinality does not exceed it.""" 

Prob={ } 

if max_losses == None or max_losses > len(vulnerable_edges): 

max_losses = len(vulnerable_edges) 
for i in range(max_losses+l): 

for x_tilde in combinations(vulnerable_edges,i): 
p=l 

for e in vulnerable_edges: 
if e in x_tilde: 

p *= phat[e] 
else: 

p *= 1-phat[e] 

Prob[x_tilde]=p 
#normalize probabilities 
totprob = math.fsum(Prob.values 0) 
for x_tilde in Prob: 

Prob[x_tilde] /= totprob 
return Prob 


def generate_xhat(vulnerable_edges, phat, xhat): 
for edges in vulnerable_edges: 
temp = random.random() 
if temp <= phat[edges]: 

xhat[edges]=1 
else: 

xhat[edges]=0 
return xhat 

def generate_xhat_from_x_tilde(vulnerable_edges,x_tilde,xhat): 
for edge in vulnerable_edges: 
if edge in x_tilde: 

xhat[edge]=1 
else: 

xhat[edge]=0 
return xhat 
#End phat_code.py 
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#Begin combo_model.py 

from _future_ import division 

import sys 

timport coopr and Solver 

from coopr.pyomo import * 

from coopr.opt import SolverFactory 

tMaster cuts 

def MasterStochCut(model, iter_num) : 

nodeterm = sum(model.b[i]*model.pi[i,iter_num] for i in model.Nodes) 
arcterm = sum(model.u[i,j]*model.mu[i,j,d,iter_num]*model.W[i,j,d] for 
(i,j,d) in model.Vuln*model.Defenses) 

return model.Z_Stoch >= nodeterm + arcterm 

def MasterDetCut(model,iter_num): 
edge_sum = ( 

sum((model.Y[i,j,'O',iter_num]*(model.c[i,j]+(model.xhat[i,j,iter_num]*model.q[ 
i,j])))for (i,j) in model.Vuln) 

+sum((model.Y[j,i,'O',iter_num]*(model.c[j,i]+(model.xhat[i,j,iter_num]*model.q 
[j,i])))for (i,j) in model.Vuln) 

+sum((model.Y[i,j,'1',iter_num]*(model.c[i,j])) for (i,j) in 

model.Arcs) 

) 

node_sum = sum(model.v[j]*model.S[j, iter_num] for j in model.Nodes) 
return model.Z_Det >= edge_sum + node_sum 

def MasterDetFlowBalance(model,node,iter_num): 

amountin = sum (model.Y[i,j,defense,iter_num] for (i,j,defense) in 
model.Arcs*model.Defenses if j==node) 

amountOut = sum (model.Y[i,j,defense,iter_num] for (i,j,defense) in 
model.Arcs*model.Defenses if i==node) 
supply = model.b[node] 

return amountOut - amountin - model.S[node,iter_num] <= supply 

def MasterDetEdgeCapacity(model,i,j,defense,iter_num): 

return model.Y[i,j,defense,iter_num] + model.Y[j,i,defense,iter_num] <= 
model.u[i,j]*model.W[i,j,defense] 

tChoose one defense for each vulnerable edge 
def MasterOneDefense(model,i,j): 

return sum(model.W[i,j,d] for d in model.Defenses) == 1 

#Limit number of defenses 
def MasterDefenseLimit(model): 

return sum(model.W[i,j,'1'] for (i,j) in model.Vuln) <= model.max_defenses 

tMaster Objective 

def MasterObjective(model): 

return model.det_wt*model.Z_Det + model.stoch_wt*model.Z_Stoch 

tEnsure name of networkx graph is Grph to run effectively 
# need networkx graph named Grph, dictionary for xhats 
def BuildComboMaster(Grph, vulnerable_edges, defenses, max_defenses, 
iterations, det_wt, stoch_wt, xhats, nodepi, edgemu, supply, capacity): 
tDefine model type 
model = AbstractModel() 

tPut model into Pyomo 

nodes = Grph.nodes() 

model.Nodes = Set(initialize=nodes) 


40 



model.Vuln = Set(within=model.Nodes*model.Nodes, 
initialize=vulnerable_edges) 
arcs = Grph.edges () 

model.Arcs = Set(within=model.Nodes*model.Nodes,initialize=arcs) 

iters = range(iterations) 

model.Iters = Set(initialize=iters) 

model.Defenses = Set(initialize=defenses) 

model.max_defenses = Param(initialize=max_defenses) 

#Node dual parameters 

model.pi = Param(model.Nodes*model.Iters, initialize=nodepi) 

#Arc capacity dual parameters 

model.mu = Param(model.Arcs*model.Defenses*model.Iters, initialize=edgemu) 

model.xhat = Param(model.Vuln*model.Iters, initialize=xhats) 

fare costs 
values= { } 

for item in Grph.edges(): 

values[item] = Grph.edge[item[0]][item[l]]['cost'] 
model.c = Param(model.Nodes,model.Nodes, initialize=values) 

#supplies 

model.b = Param(model.Nodes,initialize=supply) 

#capacities 

model.u = Param(model.Vuln, initialize=capacity) 

#PenaltY model 

model.q = Param(model.Arcs, initialize=100) 

#PenaltY unsupplied nodes 

model.V = Param(model.Nodes, initialize=10) 

model.det_wt = Param(initialize=det_wt) 
model.stoch_wt = Param(initialize=stoch_wt) 

#define objective variables 

model.Z_Stoch = Var(domain=NonNegativeReals) 
model.Z_Det = Var(domain=NonNegativeReals) 

#define design variables 

model.W = Var(model.Vuln*model.Defenses, domain=BinarY) 

#define response variables for deterministic cuts 
model.Y = Var(model.Arcs*model.Defenses*model.Iters, 
domain=NonNegativeReals) 

#Define fuel shortfall at Node 

model.S = Var(model.Nodes*model.Iters, domain=NonNegativeReals) 

#Master cut constraints 

model.MasterStochCut = Constraint(model.Iters, rule=MasterStochCut) 
model.MasterDetCut = Constraint(model.Iters, rule=MasterDetCut) 

#Det flow constraints for each iter 

model.MasterDetFlowBalance = Constraint(model.Nodes*model.Iters, 
rule=MasterDetFlowBalance) 
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model.MasterDetEdgeCapacity = 

Constraint(model.Vuln*model.Detenses*model.Iters, rule=MasterDetEdgeCapacitY) 
#Master arc defense constraints 

model.MasterOneDefense = Constraint(model.Vuln, rule=MasterOneDefense) 
#Master limit on number of defenses 

model.MasterDefenseLimit = Constraint(rule=MasterDefenseLimit) 

#Objective function 

model.MasterObjective = Objective(rule=MasterObjective, sense=minimize) 
instance=model.create () 
return model, instance 

tCreate the model instance and Solve returning values for 

#min flow path in variable called shortest 

#need list of double_nodes to create key and edge labels 

def SolveModel(instance): 

opt = SolverFactory("cplex") 

#opt = SolverFactory("glpk") 

#ask for duals and reduced costs from solver 
results = opt.solve(instance, suffixes=['dual','rc']) 
return results 

def BuildResults(instance,results): 
instance.load(resuits) 

#figure arcs defended 

defended = [] 

for key in instance.W: 

if instance.W[key].value > 0: 
defended.append(key) 

objective_value = instance.MasterObjective() 
return defended, objective_value 
#End combo_model.py 


42 



#Begin pyo_model.py 

from _future_ import division 

import sys 

timport coopr and Solver 

from coopr.pyomo import * 

from coopr.opt import SolverFactory 

tNetwork balance of flow constraint 
def FlowBalance(model, node) : 

amountin = sum (model.Y[i,j,defense] for (i,j,defense) in 
model.Arcs*model.Defenses if j==node) 

amountOut = sum (model.Y[i,j,defense] for (i,j,defense) in 
model.Arcs*model.Defenses if i==node) 
supply = model.b[node] 

return amountOut - amountin - model.S[node] <= supply 

tDefine objective equation 
def FlowCost(model): 
edge_sum = ( 

sum((model.Y[i,j,'0']*(model.c[i,j]+(model.xhat[i,j]*model.q[i,j])))for (i,j) 
in model.Vuln) 

+sum((model.Y[j,i,'0']*(model.c[j,i]+(model.xhat[i,j]*model.q[j,i])))for (i,j) 
in model.Vuln) 

+sum((model.Y[i,j1']*(model.c[i,j])) for (i,j) in model.Arcs) 

) 

node_sum = sum(model.v[j]*model.S[j] for j in model.Nodes) 
return edge_sum + node_sum 

def EdgeCapacity(model, i, j, defense): 

return model.Y[i,j,defense] + model.Y[j,i,defense] <= 
model.u[i,j]*model.w_hat[i,j,defense] 

tEnsure name of networkx graph is Grph to run effectively 
# need networkx graph named Grph, dictionary for xhats 

def build_model(Grph, xhat, w_hat, vulnerable_edges, defenses, supply, 
capacity): 

#Define model type 
model = AbstractModel() 

#Put model into Pyomo 
nodes = Grph.nodes]) 

model.Nodes = Set(initialize= nodes) 
edges = Grph.edges]) 

model.Arcs = Set(within=model.Nodes*model.Nodes,initialize = edges) 

model.Vuln = Set(within=model.Nodes*model.Nodes,initialize = 
vulnerable_edges) 

model.Defenses = Set(initialize=defenses) 

#create suffixes for dual information 

model.dual = Suffix(direction=Suffix.IMPORT_EXPORT) 

#Create model for costs 

values= {} 

for item in edges: 

values[item] = Grph.edge[item[0]][item[l]]['cost'] 
model.c = Param(model.Nodes,model.Nodes, initialize=values) 
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#Supplies 

model.b = Param(model.Nodes,initialize=supply) 

#PenaltY model 

model.q = Param(model.Arcs, initialize=100) 

#PenaltY unsupplied nodes 

model.V = Param(model.Nodes, initialize=10) 

#Capacities on vulnerable arcs 

model.u = Param(model.Vuln, initialize=capacitY) 

#Define attack variable 

model.xhat = Param(model.Vuln,initialize=xhat) 

#Define attack variable 

model.w_hat = Param(model.Vuln*model.Defenses,initialize=w_hat) 
#Define Variable we want to change 

model.Y = Var(model.Arcs*model.Defenses, domain=NonNegativeReals) 
#Define fuel shortfall at Node 

model.S = Var(model.Nodes, domain=NonNegativeReals) 

#Constraint for balance 

model.FlowBalance = Constraint(model.Nodes, rule=FlowBalance) 

model.EdgeCapacity = Constraint(model.Vuln*model.Defenses, 
rule=EdgeCapacity) 

#Objective function 

model.ElowCost = Objective(rule=ElowCost, sense=minimize) 
instance=model.create () 

return model, instance 

tCreate the model instance and Solve returning values for 

#min flow path in variable called shortest 

#need list of double_nodes to create key and edge labels 

def solve_model(instance): 

opt = SolverEactory("cplex") 

#opt = SolverEactory("glpk") 

#ask for duals and reduced costs from solver 
results = opt.solve(instance, suffixes=['dual' , 'rc']) 
return results 

def build_results(instance,results): 
instance.load(resuits) 

#figure arcs used 
shortest = [] 
edge_values = {} 
for key in instance.Y: 

if instance.Y[key].value > 0: 

shortest.append((key[0],key[1])) 

edge_values[(key[0],key[1])] = instance.Y[key].value 

objective_value = instance.ElowCost(instance) 
return shortest, edge_values, objective_value 

#End pYO_model.pY 
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